First order phase transition in the anisotropic quantum orbital compass model 
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We investigate the anisotropic quantum orbital compass model on an infinite square lattice by 
means of the infinite projected entangled-pair state algorithm. For varying values of the Jx and 
Jz coupling constants of the model, we approximate the ground state and evaluate quantities such 
as its expected energy and local order parameters. We also compute adiabatic time evolutions of 
the ground state, and show that several ground states with different local properties coexist at 
Jx = Jz- All our calculations are fully consistent with a first order quantum phase transition at this 
point, thus corroborating previous numerical evidence. Our results also suggest that tensor network 
algorithms are particularly fitted to characterize first order quantum phase transitions. 
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Introduction.- When quantum many-body systems are 
cooled down close to zero temperature, important collec- 
tive phenomena may occur [1]. A good example is pro- 
vided by transition-metal oxides, whose physical proper- 
ties have become of increasing interest in the last few 
years [2]. In these compounds the orbital degrees of 
freedom of the atomic electrons play a key role in de- 
termining properties such as metal- insulator transitions, 
high-temperature superconductivity and colossal magnc- 
toresistance. 

The paradigmatic approach to these systems is based 
on the so-called orbital compass models [3, 4], which have 
been the subject of many studies in the past both in the 
classical and quantum regimes. For these systems, Jahn- 
Tcllcr effects produce an anisotropy of the pseudospin 
couplings which is intertwined with the orientation of the 
interaction bonds. The properties of these systems have 
attracted considerable attention since they are endowed 
with symmetries that effectively reduce the dimension- 
ality of the system (the so-called dimensional reduction) 
[5, 6]. Despite of their apparent simplicity, orbital com- 
pass models are relevant in a variety of contexts, such 
as in determining the physics of Mott insulators with or- 
bital degrees of freedom [.3] and the implementation of 
protected qubits for quantum computation in Josephson 
junction arrays [7]. These systems are also candidates to 
exhibit topological quantum order [8]. Furthermore, it 
was recently shown how to simulate these models using 
polar molecules in optical lattices and systems of trapped 
ions with state-of-the-art technology [9, 10]. 

Generally speaking, the symmetries in these systems 
involve large degeneracies in their energy spectra, which 
make their numerical simulation difficult [11]. This 
fact, together with the lack of exact solutions, makes 
it hard to elucidate their phase diagrams. In this pa- 
per we use a tensor product state (TPS) [12, 13] or pro- 
jected entangled-pair state (PEPS) [14] to study the two- 
dimensional anisotropic quantum orbital compass model 
(AQOCM) and, in particular, to investigate whether its 
phase transition is of first order [15, 16] or second order 
[17]. More specifically, we use the infinite PEPS (iPEPS) 



algorithm of Ref. [18] to study the model directly in the 
thermodynamic limit. Our results provide abundant evi- 
dence in favor of a first order phase transition. 

The model.- The 2D AQOCM describes a system of 
spins 1/2 interacting on a square lattice with anisotropic 
two-body interactions as defined by the Hamiltonian 



(id) 



J: 



(1) 



is the Pauli X (Z) operator at site 
and Jx (Jz) is the coupling in the x 



where Xl'^l {Z^''^^ 
[i,j] of the lattice, 
(z) direction. 

For this model, Nussinov and Fradkin [1!)] proved that 
its Hamiltonian is dual to a plaquette model proposed by 
Xu and Moore to describe p-\- ip superconducting arrays 
such as Sr2Ru04 [17]. The influence of impurities [20] 
and of diluted lattices [21] in the model has also been 
investigated. In addition, finite temperature properties 
have been studied both in the quantum and classical ver- 
sions of the model [11, 22], and in both cases the existence 
of a low temperature ordered phase with a thermal tran- 
sition lying in the 2D Ising universality class has been 
shown. Finally, in Ref. [23] a ID version of the model 
was shown to undergo a first order phase transition. 

The Hamiltonian from Eq. (1) has also some signifi- 
cant properties in the context of quantum computation. 
For instance, the model was proven to be dual to the 2D 
cluster state Hamiltonian embedded in a magnetic field 
[24] . It was also shown to be related to certain classes of 
quantum error correcting codes where the system is used 
to codify a qubit that is robust against external local 
noise [2")]. 

Before proceeding any further, let us sketch some of 
the basic symmetry properties of the Hamiltonian H in 
Eq. (1) (see e.g. Refs [19, 25] for detailed discussions). 
Define the operators 
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where Pi acts on column i of the 2D lattice and Qj acts 
on row ). It is not difficult to check that these operators 



2 



(i) 



(ii) 



-0.40 
-0.42 
~-0.44 
-0.46 
-0.48 
-0.50| 

^-0.4018 

j -0.4019 
I 

i -0.4020 
) 

-0.4021 



0.5 



1 




D=5 D=6 

\ = 2 \ \ 

D=3 0=4 



Finite Size 
Anaiysis 



0.4 



1/D 



0.2 



0.005 

1/N 



-0.4018 
-0.4019 
-0.4020 
-0.4021 



FIG. 1: (color online) (i) Energy per link e in the AQOCM 
on an infinite square lattice obtained by using the iPEPS al- 
gorithm with D = 2,3 (results for — 4,5,6 are very similar 
to those for D = 3). The energy e has a sharp peak with 
discontinuous derivative at the phase transition. Dotted lines 
correspond to the results from Ref. [I 'l] up to 16x16 lat- 
tices using exact diagonalization and Green's function Mon- 
tecarlo (plotted with permission). Lines linking numerical 
points are a guide to the eye. (ii) Comparison at s = 0.5 of 
the energy per bond computed with the iPEPS algorithm for 
D = 2, . . . ,6 and the finite-size analysis from Ref. [lo]. Lines 
linking numerical points are a guide to the eye. 



commute with H for all the values of i and j. Impor- 
tantly, [Pi, Qj] 7^ for any i, j, and therefore operators 
Pi and Qj represent incompatible symmetries of H. Fur- 
thermore, notice that [Pi,Pii] = \fi,i' and similarly 
for Qj, and that any tensor product of operators corre- 
sponding to different columns (or rows) commutes with 
H as well. All these symmetries imply that, in the case 
of a system defined on SlJi L x L square lattice, every 
eigenstate of the Hamiltonian is at least of order 0{2^) 
degenerate. Also, whenever Jx = Jz the system is invari- 
ant under the reflection symmetry X Z , indicating 
the self-duality of the model at equal couplings [.^')]. 

The above self-duality indicates a possible phase tran- 
sition in the system at Jx — Jz- There have been sev- 
eral attempts to determine the existence and order of 
this phase transition. On the one hand, Xu and Moore 
pointed towards a possible second order quantum phase 
transition [17]. On the other hand, some approximate 
calculations seem to favor a first order transition [1.'3, IG]. 
The nature of this phase transition is, therefore, not to- 
tally understood yet. 

The method. - In this paper we use the iPEPS algorithm 
[18] to compute the ground state as well as adiabatic time 
evolutions for the AQOCM on an infinite 2D square lat- 
tice. As explained in Refs. [14, 18], the accuracy of the 
results relies on a refinement parameter that we shall re- 
fer to as D. This parameter is related to the maximum 
entanglement content that can be handled by the simu- 
lations [29]. In practice, increasing the value of D leads 
to better descriptions of the ground state, and therefore 
to more accurate estimations of the different observable 
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FIG. 2: (color online) Expected energy per lattice link for the 
ground state |*I'gs(s)) and the adiabatically evolved states 
|L(s)) and |JfI(s)), as computed with the iPEPS algorithm 
with D = 2. 



quantities. In our calculations we consider D = 2, . . . , 6 
and, without loss of generality, Jx,Jz > [■>()]. The cou- 
pling strengths in Eq. (1) can be restricted to the range 
Jx, Jz G [0, 1] and written in terms of a variable s S [0, 1] 
as Jx — cos (s7r/2) and J^ ~ sin (s7r/2). 

Let us discuss the impact that the symmetries of the 
system have in our simulations. As explained above, the 
symmetries of the AQOCM imply an infinite degener- 
acy of its ground state in the thermodynamic limit [19]. 
For instance, different ground state wave functions can 
be labelled according to the different eigenvalues of op- 
erators Pi in Eq. (2). This sort of degeneracy, however, 
does not play a significant role in our simulations since 
our representation of the quantum state by means of an 
iPEPS is, by construction, invariant under translations 
in the x and z directions [Is]. Still, our implementation 
of the algorithm could be sensitive to the two-fold degen- 
eracy caused by a simultaneous flip of all the spins. In 
practice, however, we observe that this does not happen. 
The simulations spontaneously choose either a positive 
or negative value of (Xl^'-'l) (or {Z^^'^^)) for all sites [i,j] 
away from the phase transition point [31]. 

Simulation results.- Our calculations are of two types. 
First, we have computed the ground state wave function 
^Gs(s)) of the system as a function of s and evaluated 
observable quantities on it such as energy and local or- 
der parameters. Second, we have simulated adiabatic 
time evolutions starting from the computed ground state 
I'^Gsisini)) for a given initial parameter Sini, and adia- 
batically increasing or decreasing s in the Hamiltonian 
well beyond crossing the point s = 1/2 {Jx — Jz)- These 
evolutions define two families of states, the left \L{s)) for 
Sini < 1/2, and the right |i?(s)) for Sini > 1/2. 

The ground state energy per lattice link 

e{s) ^ |i(M/G5(s)|X[^^^lx['+i'^llvI/G5(s)) 

+ ^{^SJGs{s)\Z^'''^Z^^+'-^^^\^Gsis)), (3) 

(independent of i and j) is displayed in Fig. (1). Our 
results show the presence of a sharp peak at s ~ 1/2, 
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which is compatible with the existence of a first order 
phase transition at this point. The energy per Unk in 
the adiabatically evolved states \L{s)) and \R{s)) is also 
plotted in Fig. (2). There we can see that the energy 
of e.g. \L{s)) follows the ground state energy up to the 
transition point s = 1/2. More generally, we find that, up 
to numerical accuracy, the PEPS for \L{s)) is the same 
as that for the ground state for s < 1/2 (and similarly 
for |i?(s)) in the regime s > 1/2). Therefore, 



I*GS(S)) 



\L{s)) if s < 1/2 
\R{s)) if s > 1/2 



From Fig. (2) we can also infer that state \L{s)) no 
longer corresponds to the ground state of the system for 
s > 1/2, but rather to some higher-energy excitation 
(and similarly for |-R(s)) for s < 1/2). The simulations 
of \L{s)) and \R{s)) are robust against modifying the 
rate of change of the Hamiltonian during the adiabatic 
evolution, indicating the presence of an energy gap to the 
reachable excitations. At the phase transition point both 
states |i(l/2)) and |i?(l/2)) have the same energy as the 
actual ground state |4'gs(1/2)), indicating the presence 
of two possible ground states of the system at this point. 

Importantly, these two ground states at s = 1/2 can be 
shown to be locally different, for instance by computing 
the Ising-likc order parameters 

m,(.s) = (*gs(s)|X['''^1|*gs(s)), 
m,{s) = (vI/gs(s)|^[*'^'1|*gs(s)), 



(4) 
(5) 



which arc independent of due to translation in- 

variancc. Fig. (3) shows rrix and ruz as a function 
of s, together with analogous expected values m^(s), 
m^{s), m^{s) and m^{s) for the evolved states \L{s)) 
and |i?(s)). We find that rux and rriz are both discon- 
tinuous at s = 1/2. However, such discontinuity could 
originate in a lack of resolution in s. That is, perhaps by 
considering more points around s = 1/2, the discontinu- 
ity in the order parameters would disappear, indicating a 
continuous phase transition. This possibility can be ruled 
out by noticing that e.g. m^(s) does not vanish to the 
right of the transition point (similarly, mf{s) does not 
vanish to the left of the transition point). That is, the 
two families of states \L{s)) and \R{s)), which coincide 
with the ground state to the left (respectively right) of 
s = 1/2, remain locally different at the transition point, 
where both represent possible ground states of the sys- 
tem. We interpret this fact as conclusive evidence of the 
existence in the 2D AQOCM of a first order phase tran- 
sition between the two phases characterized by vanishing 
and non-vanishing values of the local order parameters 
nix and m^. 

Let us now discuss the role played by the symmetries in 
this phase transition. Our numerical calculations using 
tensor networks have also shown that the ground states 
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FIG. 3: (color online) Expected values of the local order pa- 
rameter operators x'*'-'" and z'''-'' (in absolute value) for the 
ground state |^'gs(s)) {rrix and rriz) and the adiabatically 
evolved states \L{s)) (m^ and m^) and |i?(s)) (m^ and mf), 
obtained by using the iPEPS algorithm with D — 2,3 (results 
for D = 4,5 are very similar to those for D = 3). The lines 
correspond to the results from Ref. [16] using mean field the- 
ory after fermionization of the Hamiltonian (plotted with the 
author's permission). 



I'I'gs(s)) satisfy the eigenvalue relations Pi|^'Gs(s)} = 
|*Gs(s)) if s < 1/2 and Qj|*gs(s)) = |*gs(s)) if 
s > 1/2, regardless of the values of i and j. Thus, we see 
that the system chooses to preserve a different symmetry 
at each side of the phase transition point, namely, the Pi 
symmetry for s < 1/2 and the Qj symmetry for s > 1/2. 
Quite naturally, the system chooses to break the symme- 
try which minimizes the amount of entanglement in the 
broken ground state, while leaving the remaining symme- 
try intact. In turn, this also implies that the adiabatically 
evolved states \L{s)) and |i?(s)) are, respectively, eigen- 
states of operators Pi and Qj with eigenvalue 1 for any 
value of s. This follows from the fact that the symme- 
try of the initial state is preserved all along the adiabatic 
time evolution since the symmetry operators commute 
with the Hamiltonian for any value of s. Therefore, the 
two possible ground states at the phase transition point 
L(l/2)) and |i?(l/2)) obtained by adiabatic evolution 
preserve the Pi and Qj symmetries respectively. 

In addition, we observe that the two families of adi- 
abatically evolved states are related to each other by a 
non-local transformation, namely the duality transfor- 
mation of the model that switches the values of Jx and 
Jz in Eq. (1). More precisely, for all the computed 
values of s, these are related by a rotation \L{s)) = 
W{TT/2)w{7r/2)\R{s)}, where W^(7r/2) rotates the spin de- 
grees of freedom by an angle 7r/2 around the y-axis and 
w{tt/2) rotates the square lattice by tt/2. That it takes a 
highly non-local transformation to map \L{s)) and \R{s)) 
into each other is, again, consistent with a first order 
transition, where the two coexisting ground states are 
not expected to be connected by local perturbations. 

Furthermore, we have also computed the ground state 
fidelity-per-site diagram [26, 27, 28] for this system (not 
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shown) and have obtained results that agree with the 
typical behavior expected of a first order transition (see 
Ref. [28]). 

All the above results are compatible with those ob- 
tained using other numerical approaches. As a first 
check, we have verified that our simulations reproduce 
the results of simple series expansion calculations that 
we performed far away from s = 1/2. As can be seen 
in Fig. (1), the present results for the energy per bond 
e, computed directly for an infinite system, agree in the 
first 4 significant digits with the value obtained through 
a rough extrapolation, to the thermodynamic limit, of 
exact diagonalization and Green's function Montecarlo 
results for finite systems presented in Ref. [15]. More- 
over, as shown in Fig. (3), close to the phase transition 
point the present results for the order parameters nix 
and niz are comparable to those obtained in Ref. [l(i] 
with mean field theory after fermionization of the Hamil- 
tonian. The small disagreement, of the order of 1.5 %, 
increases with growing values of D (that is, as our results 
become more precise), which suggests that the iPEPS re- 
sults for D = 2 are already better than those obtained 
by combining fermionization with mean field theory. We 
stress that our simulations show fast convergence of the 
computed observables with the refinement parameter D 
(see e.g. Fig. (l.(n))). 

Conclusions.- In this paper we have provided fresh ev- 
idence that, contrary to what had been suggested in Ref. 
[17], the phase transition in the AQOCM on a square 
lattice is of first order. Unlike previous approaches to 
this problem, we have employed an algorithm based on 
a TPS or PEPS for an infinite 2D lattice to numerically 
compute the ground state and, for the first time for an 
infinite 2D system, its adiabatic time evolution. We be- 
lieve that our results, together with those in Ref. [15, IG], 
conclusively support the existence of a first order phase 
transition. 
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